OUDLEV’; >. LIBRARY 
NAVAL f HADUATE SCHOOL 
MONTEfvL CA 93943-5101 



Approved for public release; distribution is unlimited 



Effects of Observer Dynamics 
on Motion Stability of Autonomous Vehicles 

by 



Bulent picay 

Lieutenant J.G., Turkish Navy 
B. S., Turkish Naval Academy, 1988 



Submitted in partial fulfillment of the 
requirements for the degree of 

MASTER OF SCIENCE IN MECHANICAL ENCINEERINC 

from the 

NAVAL POSTGRADUATE SCHOOL 
June 1994 



/ 



REPORT DOCUMENTATION PAGE 



Form Approved 
0MB No. 0704-0188 



PuWic reporting burden for this collection of information is estimated to average l hour per response, including the time for reviewing instruaions. searching e«isting data sources, 
gathering and maintaining the data needed, and completing ar>d reviewing the collection of information Send comments regarding this burden estimate or any other aspect of this 
colleaion of information, including suggestions for reducing this burden to Washington Headquarters Services. Directorate tor Information Operatiorts and Reports. 1215 Jefferson 
Davis Highway. Suite 1204. Arlington, VA 22202-4302. ar>d to the Office of Management and Budget. Paperwork Reduaion Projea (0704-0188). Washington. DC 20503. 



1. AGENCY USE ONLY (Leave blank) 


2. REPORT DATE 




June 1994 



3. REPORT TYPE AND DATES COVERED 

Master’s Thesis 



4. TITLE AND SUBTITLE 

EFFECTS OF OBSERVER DYNAMICS ON MOTION STABILITY 
OF AUTONOMOUS VEHICLES 



6. AUTHOR(S) 

Okay, Biilent 



5. FUNDING NUMBERS 



7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 



PERFORMING ORGANIZATION 
REPORT NUMBER 



Naval Postgraduate School 
Monterey, CA 93943-5000 



9. SPONSORING /MONITORING AGENCY NAME(S) AND ADDRESS(ES) 



10. SPONSORING /MONITORING 
AGENCY REPORT NUMBER 



11. SUPPLEMENTARY NOTES 

The views expressed in this thesis are those of the author and do not reflect the 
official policy or position of the Department of Defense or the United States Government. 



12a. DISTRIBUTION /AVAILABILITY STATEMENT 

Approved for public release; distribution unlimited. 



12b. DISTRIBUTION CODE 



13. ABSTRACT (Maximum 200 words) 

The problem of loss of stability of marine vehicles under cross track error control in the presence of 
mathematical versus actual system mismatch is analyzed. For demonstration purposes, variations in the 
heading angle control gain are studied. Particular emphasis is placed on analyzing the effects of observer 
design on system response after initial loss of stability of straight line motion. It is shown that the dynamics 
of the observer may have a significant effect on the computed gain margin of the control system depending 
on the particular basis used. 



14. SUBJECT TERMS 

Hopf bifurcation, supercritical, center manifold theorem 
limit cycle, periodic solutions 



15. NUMBER OF PAGES 

68 



16. PRICE CODE 



17 . 



SECURITY CLASSIFICATION 
OF REPORT 



TTNrT.ASqTFTF.n 



18. 



SECURITY CLASSIFICATION 
OF THIS PAGE 



UNCLASSIFIED 



19. 



SECURITY CLASSIFICATION 
OF ABSTRACT 



20. LIMITATION OF ABSTRACT 



UNCLASSIFIED 



_LIL 



NSN 7540-01-280-5500 



Standard Form 298 (Rev 

Prevrnbrd by ANSI Std Z39-10 
298-102 



2-89) 



ABSTRACT 



The problem of loss of stability of marine vehicles under cross track error 
control in the presence of mathematical versus actual system mismatch is analyzed. 
For demonstration purposes, variations in the heading angle control gain are studied. 
Particular emphasis is placed on analyzing the effects of observer design on system 
response after initial loss of stability of straight hne motion. It is shown that the 
dynamics of the observer may have a significant effect on the computed gain margin 
of the control system depending on the particular basis used. 



Ill 



TABLE OF CONTENTS 

I. INTRODUCTION 1 

II. PROBLEM FORMULATION 3 

A. EQUATIONS OF MOTION 3 

B. COMPENSATOR DESIGN 4 

C. CALCULATION OF CONTROL GAINS 7 

D. CALCULATION OF OBSERVER GAINS 8 

E. CHARACTERISTICS OF ITAE CRITERIA 9 

III. HOPE BIFURCATION 11 

A. INTRODUCTION 11 

B. THIRD ORDER EXPANSIONS OF THE SYSTEM EQUATIONS 12 

1. Perturbation in 12 

2. Integral Averaging 14 

C. RESULTS 17 

IV. COMPENSATOR IN A DIFFERENT BASIS 19 

A. CRITICAL VALUE OF C 19 

B. THIRD ORDER EXPANSIONS OF THE SYSTEM EQUATIONS 19 

1. Perturbation in Ki^ 19 

2. Integral Averaging 23 

C. RESULTS 25 

V. CONCLUSIONS AND RECOMMENDATIONS 29 

A. CONCLUSIONS 29 

D. RECOMMENDATIONS 30 

APPENDIX A: HOPF BIFURCATION PROGRAM FOR [X,X] BASIS . . 31 

APPENDIX B: CRITICAL VALUE OF (7 FOR [Y,Y] BASIS 41 

APPENDIX C: HOPF BIFURCATION PROGRAM FOR [X,X] BASIS ... 47 

REFERENCES 57 

INITIAL DISTRIBUTION LIST 59 

iv 



DUDLEY KNOX LIBRARY 
NAVAL POSTGRADUATE SCHOOL 
MONTEREY CA 93943-5101 



LIST OF FIGURES 

2.1 Saturation in <5 5 

2.2 Step response of ITAE 9 

4.1 Ccrit versus natural frequency for 20 

4.2 Kfc^ versus Wn for diflFerent observer 27 



V 



TABLE OF SYMBOLS 



a 


dummy independent variable, or yaw rate coefficient in Nomoto’s model 


A 


linearized system matrix 


b 


rudder angle in Nomoto’s model 


c 


parameter for variance of gain and hydrodynamic coefficients 


Cent 


bifurcation value of c 


Iz 


vehicle mass moment of inertia 


K 


cubic stability coefficient 


K^, Kr, Ky 


controller gains 


£\f/, £ri £y 


observer gains 


m 


vehicle mass 


N 


yaw moment 


PAH 


Poincare- Andronov- Hopf Bifurcation 


r 


yaw rate 


R 


polar coordinate of transformed reduced system 


T 


matrix of eigenvectors of A, or limit cycle period 


V 


sway velocity 


X 


state variables vector 


Xg 


body fixed coordinate of vehicle center of gravity 


y 


deviation of the commanded path 



VI 



V 



sway force 



z stable variables vector in canonical form 

zi, Z 2 critical variables of z 

ao, Oi, 0!2 coefficients of desired characteristic equation 
real part of critical pair of eigenvalues 
j3' derivative of ^ with respect to c evaluated at Ca-u 

70) 7ii 72 coefficients of desired characteristic equation 
6 rudder angle control law 

6q linearized rudder angle control law 

£ critical difference c — Cait 

9 polar coordinate of transformed reduced system 

ip vehicle heading angle 

u imaginary part of critical pair of eigenvalues 

u}' derivative of u with respect to c evaluated at Ccrit 

u>n natural frequency 

Uno observer natural frequency 



VII 



ACKNOWLEDGEMENT 



I wish to thank my thesis advisor, Professor Fotis A. Papouhas, for his guidance 



and encouragement in this research. 



VllI 



I. INTRODUCTION 



Accurate path control of surface ships and underwater vehicles along prescribed 
geographic 2 d paths is a basic problem that is becoming increasingly important, par- 
ticularly as the missions of ocean vehicles become more complicated with strict 
requirements for performance. In order for a control law to be able to perform its 
mission in a realistic operational scenario it has to be robust enough so that it can 
maintain stability and accuracy of operations in the presence of modeling errors and 
environmental uncertainties. The robustness properties of the design are particu- 
larly important due to the unpredictable nature of the ocean environment and the 
changes in the hydrodynamic characteristics of the vehicle during turning, changes 
in the forward speed, or operations in proximity to other objects in the area. For 
these reasons, there exists a need for the analysis of the robustness characteristics 
of a particular control law design and the establishment of a rational operational 
envelope based on stability and performance criteria. Previous studies [Parsons 
and Cuong (1977)] showed that gain adaption is highly desirable due to changes 
in the linearized vehicle hydrodynamics with different operation conditions, such 
as depth under keel. The resulting adaptation scheme [Parsons and Cuong (1980)] 
required significant vehicle motion, which may be undesirable when operating in 
restricted waters, or in object recognition and localization tasks. Integral control 
techniques [Parsons and Cuong (1981)] proved quite effective, but neglected the 
nonlinear behavior of the vehicle, which becomes very important at low speeds and 
hover operations. Model based compensators exhibit robust behavior under con- 
ditions of parameter uncertainty, which is as good as the classical linear quadratic 
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regulators for lineaj output feedback systems [Healey (1992)]. Alternatively, sliding 
mode controllers exhibit very robust characteristics given an estimate of the param- 
eter uncerttiinty and/or disturbances [Papoulias and Healey (1992)], [Yoerger and 
Slotine (1985)]. Sliding mode control, however, does not offer an infinitely robust 
design and it suffers from a series of bifurcation phenomena and loss of stabifity 
unless proper care is exercised [Papoulias (1991)]. 

In this work we analyze the problem of loss of stability of a marine vehicle 
under cross track error control in the presence of mathematical versus actual system 
mismatch. For demonstration purposes, variations in the heading angle control 
gain are studied. Previous studies [Oral (1993)] concentrated on system response 
assuming perfect and complete state measurement. Particular emphasis in this work 
is placed on analyzing the effects of observer design on system response after initial 
loss of stability of straight line motion. The main loss of stability cases analyzed 
here occur in the form of generic bifurcations to periodic solutions [Guckenheimer 
and Holmes (1983)]. We use center manifold reduction techniques and averaging 
in order to capture the stability properties of the resulting limit cycles [Chow and 
Mallet-Paret (1977)]. It is shown that the dynamics of the observer may have a 
significant effect on the computed gain margin of the control system depending 
on the particular basis used. All computations in this work axe conducted for the 
NPS autonomous underwater vehicle [Bahrke (1992)] and all results are presented in 
standard dimensionless quantities with respect to vehicle length, 7.3 ft, and nominal 
forward speed, 2 ft /sec. 
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II. PROBLEM FORMULATION 



A. EQUATIONS OF MOTION 

The linear maneuvering equations of motion of a marine vehicle in the hori- 
zontal plane are written in dimensionless form as, 

m(i> -H r -I- xgt) = iff- -t- Y{,v -1- -|- Ys8 (2.1) 

I^r -I- mxG(i> -I- r) = N^r N{,i/ NrV -|- N^i/ + NsS (2-2) 

where all symbols are explained in the nomenclature. Equations (2.1) and (2.2) can 
be used to derive a second order transfer function between the rudder angle 6 and 
yaw rate r. For low frequency maneuvering motions this second order equation can 
be approximated by expanding in Taylor series and keeping the first order terms 
only. The result is 

r = ar + bS (2.3) 

Equation (2.3), which is sometimes referred to as Nomoto’s first order model, is 
particularly useful in control system design since no sway velocity feedback is neces- 
sary. This equation predicts linear variation of the steady state turning rate versus 
rudder angle. In reality, the r-S curve has characteristics of softening spring mainly 
due to speed loss during turning. To account for this a modified version of equation 
(2.3) is used, 

r = ar + + bS (2-4) 

where is usually determined from steady state results. Finally, the model is 
complete by the incorporation of the kinematic equations, 

= r (2.5) 
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y = sin ^ 



( 2 . 6 ) 



where ^ is the vehicle heading, and y is the cross track error off a desired straight 
line path. 

B. COMPENSATOR DESIGN 

In control theory it is known that the eigenvalues of the controller are not 
affected by the eigenvalues of the observer. This allows us to design the controller 
and observer separately which is known as the separation principle. The combination 
is called a compensator. 

Equations (2.3), (2.5), and (2.6) govern the steering control of the model used 
in this section. The control law can be expressed «is, 

S = ^sat tanh (2'7) 

where around the nominal state ^ = r = y = 0we have 

So = K^^ + Krr->rKyy ( 2 . 8 ) 

S is the rudder angle and Kti and Ky are the control gains of the system. 
The linear control law is Sq. The rudder angle 6 is defined by a h)q)erbolic tangent 
function to include the saturation to our problem as shown in Figure 2.1. Saturation 
occurs at which is the saturation limit generally taken as 0.4 rad. 

The linearized form of equations of motions in the vicinity of^ = r = y = 0 

are. 



= r 


(2.9) 


r = ar + 6Jo 


(2.10) 


y = ^ 


(2.11) 
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i ^ 




6 . 



Figure 2.1: Saturation in S. 



These equations can be expressed in state space form as 



X = AX + Bu 



where 






r 

y J 



is the state vector, 



A = 



0 1 0 

0 a 0 

1 0 0 



0 

b 

0 



is the open loop dynamics matrix and 

B = 

is the control distribution vector. 

The observer equations are 

X = AX + Bu + L{Y- CX) 

where X is the estimate of X, Y is the output of the system Y = y, and C 
sensor vector C = [0 0 1]. 



( 2 . 12 ) 



( 2 . 13 ) 

is the 



0 



The error in the estimate of X is defined by 



k = X-k (2.14) 

Using equations (2.12), (2.13), and (2.14) we can obtaun 

k = {A- LC)X (2.15) 



We can rewrite equations (2.9), (2.10), and (2.11) in the form of 



X = 



■ ■ 




■ 0 


1 


0 ■ 




■ ^ ■ 




■ 0 ■ 


r 


= 


0 


a 


0 




r 


+ 


b 


. y . 




1 


0 


0 




. y . 




0 



and 



The observer gains are. 



r = [0 0 1] 



r 

y 



L = 



iq/ 

ir 



'■y J 



After performing the matrix operations we obtain 









■ 0 


1 


—£qi 




' ^ ■ 


= 


f 


= 


0 


a 






f 




. y . 




1 


0 


-(y \ 




. y . 



Using equation (2.13) we can rewrite equation (2.8) as follows. 
So = /C*(^ - 4') + Kr{r -f) + Ky{y - y) 



(2.16) 



(2.17) 



(2.18) 



(2.19) 



Finally, we can write our compensator equations in the form 





■ ■ 




’ 0 


1 


0 


0 


0 


0 ■ 




■ ^ • 






r 




0 


a 


0 


hcXq/ 


bKr 


bKy 




r 


' X ' 




y 




1 


0 


0 


0 


0 


0 




y 


k 




? 




0 


0 


0 


0 


1 


— £q/ 










f 




0 


0 


0 


0 


a 


-4 




f 




. y . 




, 0 


0 


0 


1 


0 


-4 J 




. y . 
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If we look at the matrix ceirefully we will see that it is in the form 



' ^ ' 




' A- BK 


BK 


k 




0 


A-LC 



which has the following characteristic equation, 

det[A -BK- s/] dei[A -LC-sI] = 0 



This indicates that the dynaimics of the observer are completely independent of the 
dynamics (eigenvalues) of the controller. Thus K and L can be designed separately. 



C. CALCULATION OF CONTROL GAINS 



A is the Jacobian matrix of the system 

0 1 



A = 



0 



bKi, a + bKr bKy 
1 0 0 



The characteristic equation of the matrix A is 



- (a + bl<r)\^ - blUX - bKy = 0 
If the desired characteristic equation has the general form 



+ ckjA^ + ckxA -|- 0(0 = 0 



the control gains can be found as 



Kif, 



Kr 

Ky 



£1 

b 

Q2 + a 
b 

QO 

b 



The desired characteristic equation can be written with respect to the desired natural 
frequency and some optimum coefficients. The ITAE criterion for a third order 
equation is 



+ 1.75tyn5* + 2.\bw\s + = 0 
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where u;„ is the desired controller natural frequency. 



Therefore the control gains can be calculated for a given natural frequency, as 

Qi = 

02 = 1.75u;„ 
ao = 



D. CALCULATION OF OBSERVER GAINS 



If we define A as the Jacobian matrix of the system 

0 1 -iit 



A = 



0 

1 



a —if 

0 -e, 



y J 



the characteristic equation of the matrix A is 

+ {iy — a) + — aiy)X + — aCif) = 0 

If the desired characteristic equation has the general form 

+ 72A^ + 7iA + 7o ~ 0 
the observer gains can be found as 



= a + 72 

= a£y + 7i 

£r = a£iH + 7o 

Applying the ITAE criteria, observer gains can be calculated for a given natural 
frequency as 



71 = 2.15u;2(, 

72 = 1.75u;„o 
7o = wlo 
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Normulizcil Amp ( y^wn^ 3 ) 



where Wno is the observer natural frequency. 



E. CHARACTERISTICS OF ITAE CRITERIA 

In the Ccilculations of gains we applied the ITAE criteria. If we look at Figure 
2.2 for the step response of ITAE, we see that the response gets faster as the natural 
frequency increases. For example, the settling time is 10 normalized seconds, or 10 
seconds for = Ij 1 second for = 10, and so on. 




Normalized Time ( t*wn ) 



Figure 2.2: Step response of ITAE. 
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III. HOPF BIFURCATION 



A. INTRODUCTION 

An important quantity in assessing the robustness of a particular control law 
design to parameter variations and unmodeled dynamics is the gain margin. This is 
defined as the extent to which changes can be inflicted on the system gain without 
loss in stability, to this end, we assume that the heading error gain is multiplied 
by a constant C. By definition, of Hopf bifurcation occurs when a pair of complex 
conjugate eigenvalues cross into the right hand half-plane. When this occurs the 
system will deviate from a steady solution in an oscillatory manner. This deviation 
can be either supercritical or subcritical [Seydel (1988)]. As the parameter C crosses 
the critical value, one pair of complex conjugate eigenvalues of the linear system 
matrix crosses transversely the imaginary axis. Locally, cis C approaches Ccnt> the 
periodic solutions are located on the two dimensional Eucledian plane spanned by the 
eigenvectors of the Jacobian matrix of the system which corresponds to the critical 
pair of eigenvalues. In this chap ' stability properties of the periodic solutions are 
established. In order to establish those properties the main nonhnear terms that 
dominate the system are isolated. Center manifold theory is used to reduce the 
flow to a two dimensional manifold. The method of averaging is then applied to the 
reduced system. 

The critical value of c for stability of straight line motion remains the same as 
[Oral (1993)], which is 

Cent = 0.2658 
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This is because the dynamics of the controller are independent from the dynamics 
of the observer as explained in Chapter II. 

B. THIRD ORDER EXPANSIONS OF THE SYSTEM EQUATIONS 
1. Perturbation in 

In the previous chapter we worked on the linear system. Now we are go- 
ing to introduce the nonlinear terms to our compensator. In this case the equations 
of motion are 



^ = r 

r = ar + + bS 

y — sin ^ 



where 



S = • tanh 



(£) 

^0 = CIU{^-'^) + Kr{r-f) + Ky[y-y) 



or in compact form, 



X = /(x) , X = [^, r, y, 'I, f, t/]' 



This system can be written in the form 



X = AX + g{x) 



(3.1) 

(3.2) 

(3.3) 

(3.4) 

(3.5) 

(3.6) 

(3.7) 



A is the Jacobian matrix of /(x) evaluated at X = 0, and g{x) contains all nonlinear 
terms of Equations (3.1), (3.2), and (3.3). Taylor expansion of the nonlinear terms 
about the equilibrium, where we keep the first non-vanishing nonlinear coefficients 
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only, gives 



sin^ 

6 



* _ 

6 




(3.8) 

(3.9) 



Substitution of Equations (3.8) and (3.9) into Equations (3.1), (3.2), and (3.3) gives 
us the A matrix in Equation (3.7) as follows. 



A = 



0 

hcIC^ 
1 
0 
0 
0 

The nonlinear parts are. 



1 

a^-hKr 

0 

0 

0 

0 



g{^) = 



bKy - 
0 
0 
0 
0 



0 

hcJCijf 

0 

0 

0 

1 



asr^ — 



6 

0 

0 

0 



'sat 



0 

-bKr 

0 

1 

a 

0 



——/fl 



0 

-bKy 

0 

— 

-£r 

-I 



y J 



(3.10) 



(3,11) 



If we introduce the transformation matrix (T) of eigenvectors of A evaluated at the 
bifurcation point. 



J, 

II 


ij = 1,2, 3, 4, 5, 6 


(3.12) 


T-^ = [no] 


i,j = 1,2, 3, 4, 5, 6 


(3.13) 


the linecur change of coordinates. 






X = Tz , 


2 = T~^x 


(3.14) 



transforms system (3.7) into its normal form 

i = T~^ATz + T-^g{Tz) (3.15) 
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With this transformation we get, 



r-MT = 



0 

two 

0 

0 

0 

0 



-two 0 
0 0 
0 
0 
0 
0 



0 
0 

Pi 0 



0 

0 

0 



P2 

0 

0 



0 

0 

0 

0 

Pz 



0 

0 

0 

0 

0 



0 P 4 I 



(3.16) 



Using center manifold theory we get, as shown in [Chow and Mallet-Paret (1977)], 



9 = 



0 

^21 •Z? + ^222^2 + ^’ IZZ \ Z 2 + 12 \ Z \ Z \ 
^ZlA. + ^ 32^2 + ^ZZZ\Z2 + t^Zizl 
0 
0 
0 



Substitution of Equation (3.14) into Equation (3.7) yields. 



zi = -tWoZ2 + rnzl + ri2ZiZ2 + ri3Ziz| + ri4z| 



(3.17) 



(3.18) 



Z 2 = tWoZi + T 2 izf + r22zlz2 + rTzZiz] + r244 (3.19) 



2. Integral Averaging 

We write Equations (3.18) and (3.19) in the form 



Zl = -tWoZ2 + Fi{zi, Z2) , 
Z2 = IWoZl + F2 (zi,Z2) 



(3.20) 

(3.21) 



If we introduce polar coordinates in the form. 



z\ = R cos 9 , Z 2 = R sin 9 



(3.22) 



Equations (3.20), (3.21) result in 



R = Fi{R,9) cos9 + F2{R,9) sin9 
R9 = two7? + ^ 2 ( 7 ?, ^)cos^ — Fi(i?, 0)sin^ 



(3.23) 

(3.24) 
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Equation (3.23) yields 



R = P{e)H^ 



(3.25) 



where P{9) is a 27r-periodic function in the angular coordinate 6. If Equation (3.25) 
is averaged over one cycle in 9, we get an equation with constant coefficients, 



R = KI^ 



(3.26) 



where, 

if = ^ jf''p(0).d0 (3.27) 

Equation (3.27) is simplified after evaluation of the integral as, 

^ ~ 8 ^’' 24 ] (3.28) 

where the coefficients are as follows. 

Til = ni2^2i + 

Ti 3 = rii2^24 + ^13^34 

T22 = ^22^23 + ^23^33 

T24 = Tl22^22 + ^23^32 

where the ni 2 , ni 3 , 7222 , and 7223 are the elements of the inverse of transformation 
matrix T. The values of the coefficients are given by the following expressions 



^21 



222^2222 



+ I<rmh + KyTnli + Sc^KlKrmhmsi 



+3c^K^Kym\iTn6i + ScKit K^mhm4i + 3cKq/KyTn\im4i + 3KrKym\^m^i 
+3/<'^/<'j,7725j77261 + 6c/<'*/C/<"ym4i 7725i7726i] 



15 



^23 — 



(■22 = o.j,m \2 — b( + A"r "^52 + ^ 1^62 + K\Krm\ 2 m ^2 

+ 3 c ^ K ^ KyTn \ 2 Tn 62 + ^ cK ^ Klm \ 2 m ^2 + " icK ^ K ^ m \ 2 m ^2 + ‘^ KrK ^ m \ 2 mf ,2 
+ ZK ^ Kyml 2 m &2 + 6 cK ^ KrKym 42 m 52 mQ 2 \ 

3a3Tnlim22 — ^<?K%m\^m42 + 3K^m\im32 + ^K^rn\^m&2 

•\- Z (? K%Kr (m 4 im 52 + 2 m 4 im 42 mnj + Z <? K%Ky (m 4 im 62 + 2 m 4 im 42 m 6 i) 
-\- ZcKif , K ^ ( Tnl ^ m 42 + 2 m 5 im 52 m 4 i) + ZcK^nKl { rn \ im 42 + 2 m 6 im 62 m 4 i) 
+ 3 ATrA'^ (m 6 im 52 + 2 m 6 im 62 m 5 i) + ZK^Ky (m 5 jm 62 + 2 m 5 im 52 m 6 i) 
+ 6 cKil,KrKy {Tn 4 imsim 62 + (014177152 + 7774277751) 7T76 i)] 

3037772177722 ~ [ 3 c^A'^ 77741 77742 + 3 A'^ 7775 i 77752 + 3 A '^ 77761 77762 

+ 3 c^A'|A'r (7774277751 + 2777417774277752) + Zc ^ K\Ky (7774277761 + 2777417774277762) 
+ 3 cA'* (7775277741 + 2777517775277742) + ZcKil/Kl (7776277741 + 2777617776277742) 

+3Ar (7776277751 + 2777617776277752) + ZK^Ky (7775277751 + 2t775i 7775 2 77762) 

+6cA»ArAj, (777427775277761 + ( 7774 i 77752 + 7774277751 ) 77762 )] 

, 1 

4 l = -T 777 



^24 — 



iz 2 = --m 



3 

12 



(33 = -^Olii 777 1 2 



6-11 
1 
6 
1 
2 
1 



(34 = -2»^iimi2 



k = 



ZSi 



sat 
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C. RESULTS 



The Vcilue of K is important for us to determine whether the bifurcation is 
supercriticcil {K < 0) or subcriticad {K > 0). In this study we wanted to see the 
effects of observer dynamics to our system, especially the value of K. To do that 
we used the Fortran code (Appendix A) for the numerical results. The result was 
that the value of K was not affected by changes in the observer natural frequency. 
The reason for this can be traced back to the definition of A", Equation (3.28). It 
can be seen that in the definitions for and iij only the first two eigenvectors m,i, 
mi 2 for t = 1,2,... 6 of the A matrix. Equation (3.10) appear. Therefore, we have 
to show that these remain the same for all observer natural frequencies. 

This A matrix. Equation (3.10), actually consists of 4 block matrices, each 
3x3, whidi are the same as shown in Equation (2.22). Let us denote the A matrix 
as follows. 



A = 



A B 

0 C 



The eigenvalues of A Ccin be computed by 



A -XI B 
0 C-XI 



= 0 



or 

\A - A/| • |C - A/| = 0 



We can group the eigenvalues in two different groups: Aj,,- for i = 1,2,3 are the 
eigenvalues of >1 an A 2 ,,- for i = 1,2,3 are the eigenvalues of C. The eigenvectors 
associated with these eigenvalues can be found as follows. 

For A = Ax^i, 



■ A - Ai.,/ 


B 


1 f 1 _ 


■ 0 ■ 


0 


1 


1 [ t^2 j 


0 
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and 



M-Ai,i/]bi| + [5][v2] = 0 

[0][vi] + [C — Ai,,/][u 2 ] = 0 

Since Ai^, is an eigenvalue of A and the eigenvalues of A and C are distinct, the 
matrix [C — Ai,,-7] is nonsingular which means that [^ 2 ] = 0. Therefore, we get 

[>l-Ai,/]H = 0 

which means that Vi is an eigenvector of A. Therefore, the first three eigenvectors of 
A axe the same as the eigenvectors of A and they are independent of the dynamics 
of the observer. Of course, the remaining three eigenvectors of A depend on the 
oberver natural frequency, but, as we pointed out earlier, none of these appear in 
the definition of the nonlinear stability coefficient K. 
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IV. COMPENSATOR IN A DIFFERENT BASIS 



A. CRITICAL VALUE OF C 

If we look at Equation (3.6), we C 2 ui see that the basis for our system was 

01^ A A 

[A'jX]. Now we are going to represent our system in [X,X] basis where X is 
the estimate of X. In this compensator basis the critical value of C in Equation 
(4.3) is no longer constant. Therefore, we used a Fortran code (Appendix B) to 
calculate the critical C values for different observer natural frequencies. A plot 
of these critical C values for different observer natural frequencies can be seen in 
Fig. 4.1. The observer natural frequencies are in nondimensional seconds whereas 
the control natural frequencies are normadized with respect to the corresponding 
observer natural frequencies. The system is unstable for values of C less than the 
critical value. 

B. THIRD ORDER EXPANSIONS OF THE SYSTEM EQUATIONS 

1. Perturbation in Kq/ 

In the previous chapters we worked on the [X, X] basis. Now we are 
going to represent our system in the new basis, which is [X,X], where X is the 
estimate of X. Our equations of motion were, 

^ = r 

r = ar + a^r^ + bS (4.1) 

y = sin ^ 
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ccRir 



CCRTT vs WN(nonnali2ed) 




Figure 4.1: Cent versus natural frequency for K^. 

where 

S = 8^ tanh (4-2) 

So = CK<fi^ Krf + Kyy (4.3) 

or in compact form, 

A- = /(i), ,Y= [t.r.y.i.r,!/]’' (4.4) 

This system can be written in the form 

.Y = AX + ( 7 ( 1 ) (4.5) 



20 



A is the Jacobian matrix of /(x) evaluated at X = 0, and ff(x) contains all non- 
linear terms of Equation (4.1). Taylor expansion of the nonlinear terms about the 
equilibrium, where we keep the first non-vanishing nonlinear coefficients only, gives 



sin ^ 
S 





(4.6) 

(4.7) 



Substitution of Equations (4.6) and (4.7) into Equation (4.1) gives us the A matrix 
in Equation (4.5) as follows 



■ 0 


1 


0 


0 


0 


0 


0 


a 


0 


hcKijf 


bKr 


bKy 


1 


0 


0 


0 


0 


0 


0 


0 




0 


1 




0 


0 


£r 


bcKifi 


bKr 


-ir + bK, 


. 0 


0 




1 


0 


-^v 



The nonlinear parts are. 



^(x) 



0 



asr^ — 



^ c3 

3e. ° 



6 

0 



3SLr 
0 



(4.8) 



(4.9) 



If we introduce the transformation matrix (T) of eigenvectors of A evaluated at the 
bifurcation point. 



T = [m.y] 


i,j — 1,2, 3, 4,5,6 


(4.10) 


T-^ = [nd 


ij = 1,2, 3, 4, 5, 6 


(4.11) 


the linecir change of coordinates. 






X = Tz , 


z = r-*x 


(4.12) 
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ticinsforms system (4.5) into its normal form 



z = T-^ATz + T-^g{Tz) 



-1 . 



(4.13) 



At the bifurcation point 



T~^AT = 



with uiQ > 0 and P, < 0. 



0 

Wo 

0 

0 

0 

0 



—Wo 0 

0 0 



0 

0 

0 

0 



0 
0 

Pi 0 



0 

0 

0 



0 0 
0 0 



P2 

0 

0 



0 

0 

Pz 



0 

0 

0 



0 Pa 



(4.14) 



Using center manifold theory we get, as shown in [Chow and Mallet- Paret 



(1977)1, 

0 

^ 212^1 + ^ 222^2 + ^ 232 ^ 1 2^2 + ^ 242 : 1^2 
_ ^ 31^1 + ^32^2 + ^33ZjZ 2 ~h t3AZ\z\ 

^ ~ 0 

^51^1 + ^52^2 "I" ^53^1 ^2 + ^54^1^2 

0 

Substitution of Equation (4.12) into Equation (4.5) yields, 

zi = -W 0 Z 2 + rnZi + ri2ZiZ2 + rizzizj -f- ri4Z2 

Z2 = Ulo^l + r2izf -h T22^?^2 + ’’23^1^2 + ’'24^2 



(4.15) 



(4.16) 

(4.17) 



where the terms r,j are evaluated by a Fortran code (Appendix C). 

For values of C close to its critical value. Equations (4.16) and (4.17) 

become. 



zi = a'ezi - (too + w'e)z2 + ruz\ -h ri2ZiZ2 -h riaZiZj -I- r^zf (4.18) 

Z2 = (too + w'e)Zi -f a'eZ2 -f r2iZi -f T22z\z2 + T2zZiz\ + T24z| (4.19) 



where e is the difference between C and its critical value C*. The terms a' and w' 
denote the derivative of the real and imaginary part, respectively, of the critical pair 
of the eigenvalues with respect to C evaluated at C* 
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2. Integral Averaging 

We write Equations (4.18) and (4.19) in the form, 

ij = a‘ezi — (u;o + w' e)z 2 + i^i( 2 i, 22 ) (4.20) 

Z 2 = (iwo + w'e)z\ + a'ez 2 + 1 ^ 2 ( 21 , 22 ) (4-21) 

If we introduce polar coordinates in the form, 

zi=Rcos6, Z 2 = Rsin6 (4.22) 

Equations (4.20), (4.21) result in 

R = a*eR + Fi{R, 6) cos 6 + F2{R,0) sin 9 (4.23) 

R9 = {wo + w'e)R+ F 2 {R, 9) cos 9 — Fi{R, 6) sin 6 (4.24) 

Equation (4.23) yields 

R = q'eR + P{9)1^ (4.25) 



where P{9) is a 27r-periodic function in the cingular coordinate 9. If Equation (4.25) 
is averaged over one cycle in 9 [Chow and Mallet-P£iret (1977)], we get an equation 
with constant coefficients, 

R = a'eR + KI^ (4.26) 

where 

Equation (4.27) is simplified after evaluation of the integral as, 

K = — [3rn + ^3 + T 22 + 3r24] (4.28) 

where the coefficients are ais follows 
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” 12^21 + ” 13^31 + ” 15^51 



rn = 

ri 3 = ^ 12^24 + »^ 13^34 + ” 15^54 

T 22 = ^ 22^23 + ^ 23^33 + ” 25^53 

’’24 = ”22^22 + ”23^32 + ”25^52 

where the ni2, ni3, nis, n22> ”23, and 7125 axe the elements of the inverse of trans- 
formation matrix T. The values of the coefficients ^21, ^22, ^23, ^24, ^31, ^32, ^33, and 
£34 are the same as in Chapter III. The values of the other £ij coefficients are given 
by the following expressions: 

£51 = —f>( [c^Klmli -I- Kfml^ -|- Kymh + 3 c^KlKrml^msi 

•\-3c^K^KymhmQi -|- ScKq,K^mlim4i + ZcK^Kym\^m4i -f- SKrK^mhmsi 
+3K^ KymliTTlQi + 6 c/<'*A'r/f'ym 4 im 5 iT 726 l] 

£52 = — 77142 + /C’”S2 + ^’^v’”s2 + 3 c^^|-^^r’” 42’^52 

-l- 3 c^A'jA'y 77 l 42’”62 + 3cKq/ K^ml2Tn42 -1- 3 cA'*A'y 772 g 2»”42 + 3 A'rA'^m 62’”52 

-f- 3 /f'^Ay 7 n 5277 l 62 + QcKif KrKym42Tn32me2\ 

£53 = —i>r [30^/^^7724177142 -|- 3 A^ 77 l 5 j 77252 -|- 3 A'y 77 lsi 77 l 62 

+ 3 c ^ I<lKr (m 4 i 7 n 52 -I- 2771417724277151) -I- 3 <^ K\Ky { rn \ imQ 2 + 2772417724277261) 

-|-3cA',t (»725im42 + 277251 7725277241) + ^ cKq / K ^ (»”6l”^42 + 277261 77262m4i) 

-|-3Ar Ay (m|i 77252 + 2 t 726 i 77262 ”251) + (’”51 ”162 + 27725i772527726l) 

+6cKq,I<rKy (772417725177262 + (>”41»”52 + »”42»”5l) »”6l)] 
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^54 = \Sc^K^m4iml2 + + S/fymeimgj 

-]-3c^K^Kr 42^51 + 2TTi4iTn42Tn$2^ "I" Z(^K^Ky (rn^2^'^i “I" 27714177142 * 7162 ) 

+ ZcKq / K ^ (t 7 i| 2^41 + 2771517715277142) + ScKij/Ky (m 62 ”l 41 + 277I6177I6277I42) 
+3KrI<l (7716277151 + 2t716i 7716277152) + ZK^Ky (t 71^2”^61 + 277l5i77l5277l62) 

+6cA'* A'rii'j, (771427715277161 + ( 77 I 4177 I 52 + 77l4277l5i) 77162)] 

C. RESULTS 

Existence and stability of limit cycles can be determined by analyzing the 
equilibrium points of the averaged Equation (4.26), which correspond to periodic 
solutions in zi, Z 2 as can be seen from Equation (4.22). From Equation (4.26) we 
can easily see that: 

1. If q' > 0, then 

(a) if K > 0, then unstable periodic solutions co-exist with the stable equi- 
librium for £ < 0, and 

(b) if K < 0, then stable periodic solutions co-exist with the unstable equi- 
librium for £ > 0. 

2. If q' < 0, then 

(a) if K > 0, then unstable periodic solutions co-exist with the stable equi- 
librium for £ > 0, and 

(b) if K < 0, then stable periodic solutions co-exist with the unstable equi- 
librium for £ < 0. 
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We refer io K < 0 as the supercritical, aind K > Q as the subcritical PAH bi- 
furcation. In the supercritical case, after the equilibrium state loses its stabihty 
the system converges to a stable periodic solution with amplitude which increases 
continuously as the difference e is increased. 

In the subcritical case, however, before the equilibrium state loses stability, its 
domain of attraction becomes very small since it is bounded by the amplitudes of the 
unstable limit cycles. In such a ccise, an initial disturbance of sufficient magnitude 
can throw the system off the nominal path even before its domain of attraction 
has completely shrunk to zero. As the nominal equilibrium becomes unstable, the 
system jumps to a different state of motion with a locally, at e = 0, discontinuous 
increase in the amplitude [Papoulias (1993)]. 

In our ccise, the value of a' is always negative, which can be seen easily from 
Figure 4.1. If we look at the nature of the curve for the critical value of C for 
different natural frequencies, we will see that as the value of critical C decreases for 
the same natural frequency, the system becomes unstable. 

After using a Fortran code (Appendix C) we observed that the nonlinear sta- 
bility coefficient K depends on observer dynamics. Figure 4.2 shows that for a given 
control design, the observer must be as responsive as possible to ensure negative K 
(stable limit cycle). On the other hand, for a given observer design, if the control 
law is too slow we get subcritical behavior (unstable limit cycle). 
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K vs WN 




Figure 4.2: versus u;„ for different observer Wn. 
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V. CONCLUSIONS AND RECOMMENDATIONS 



A. CONCLUSIONS 

An investigation of the nonlinear dynamic response characteristics of a marine 
vehicle has been presented. Particular emphasis in this work was placed on analyzing 
the effects of observer design on system response after initial loss of stabihty of 
straight line motion. Bifurcation theory techniques were utilized in order to assess 
that behavior. The main conclusions of this work can be summarized as follows. 

1. There exists a critical point for a certain combination of system gains and 
system parameters for stabihty of straight hne motion. The loss of stability 
occurs generically in the form of Poincare-Andronov-Hopf bifurcations. As 
the parameter crosses its critical value, a family of periodic orbits, self sus- 
tained oscillations develops. Center manifold reduction and integral averaging 
techniques were used in order to establish the direction of the bifurcation and 
stability of the resulting periodic solutions [Papouhas, Oral (1993)]. 

2. For [X, X] basis the critical point does not depend on the observer dynamics 
(separation principle). The nonhnear stabihty coefficient K was not influenced 
by observer dynamics. The previous reduction process shows that K depends 
on the first two eigenvectors of the 6x6 matrix A. Matrix algebra shows that 
these eigenvectors are associated only with the controller dynamics. 

3. For [X, X] basis the critical value depends on observer dynamics. For a given 
control design, the observer must be as responsive as possible to maximize the 
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region of stability. On the other hand, for a given observer design, the control 
must be as slow as possible to maximize region of stabihty. 

4. The nonlinear stability coefficient K depends on observer dynamics for this 
basis. For a given control design, the observer must be as responsive as possible 
to ensure negative K (stable limit cycles). In this benign loss of stability 
the resulting periodic solutions are continuous single-valued functions of the 
parameter distance from its critical value. On the other hand, for a given 
observer design, if the control law is too slow we get subcritical behavior 
(unstable limit cycles). In such a case, the periodic solutions develop with 
what appears to be a discontinuous increase in the amplitude of oscillations 
[Papoulias, Oral (1993)]. 

B. RECOMMENDATIONS 

The differences between the two bases with respect to robustness properties of 
the system have to be analyzed. 
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APPENDIX A 

HOPE BIFURCATION PROGRAM FOR [X,X] BASIS 



PROGRAM HTKPSI 
HOPF BIFURCATIONS 
NOMOTO'S FIRST ORDER MODEL 

CALCULATIONS FOR K AND CCRITICAL IF KPSI CHANGES W/ C 
C234567 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

C 

DOUBLE PRECISION K1 , K2 , K, LPSI, LY, LR, IZ , L, 

& MASS, NV, NR, NVDOT, NRDOT, NDRS , NDRB, KPSI , KR, KY, K3 , 

& L2 1 , L22 , L2 3 , L24 , L3 1 , L32 , L3 3 , L3 4 , LSI , L52 , L53 , L54 , 

Sc M11,M12,M13,M14,M15,M16,M21,M22,M23,M24,M25,M26, 

Sc M31 , M32 , M33 , M34 , M35 , M3 6 , M41 , M42 , M43 , M44 , M45 , M46 , 

Sc M51,M52,M53,M54,M55,M56,M61,M62,M63,M64,M65,M66, 

Sc Nil , N12 , N13 , N14 , N15 , N16 , N21, N22 , N23 , N24 , N25 , N2 6 , 

Sc N31 , N32 , N33 , N34 , N35 , N3 6 , N41, N42 , N43 , N44 , N45 , N46, 

Sc N51 , N52 , N53 , N54 , N55 , N56 , N61 , N62 , N63 , N64 , N65 , N66 

C 

DIMENSION AMAT (6,6),T(6,6), TINV (6,6), FVl ( 6 ) , IVl ( 6 ) , YYY (6,6) 
DIMENSION WR ( 6 ) ,WI(6) ,TSAVE(6,6) ,TLUD(6,6) , IVLUD ( 6 ) ,SVLUD(6) 
DIMENSION ASAVE ( 6 , 6 ) , A1 ( 6 , 6 ) , A2 ( 6 , 6 ) 

OPEN (11,FILE='AKPSI.MAT' , STATUS= 'NEW' ) 

C 

WEIGHT=435 . 0 

IZ =45.0 

L =7.3 

RHO =1.94 

G =32.2 

XG =0.0104 

MASS =WEIGHT/G 

MASS =MASS/ (0 .5*RHO*L**3) 

IZ =IZ/ (0.5*RHO*L**5) 

XG =XG/L 
YRDOT =-0.00000 
YVDOT =-0.03430 
YR =+0.00000 
YV =-0.10700 
YDRS =+0.01241 
YDRB =+0.01241 
NRDOT =-0.00047 
NVDOT =-0.00000 
NR =-0.00390 
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NV =-0.00000 
NDRS =-0.337*YDRS 
NDRB =+0.283*YDRS 
DH = (IZ-NRDOT) * (MASS-YVDOT) - 
& (MASS*XG-YRDOT) * (MASS*XG-NVDOT) 

All=( (IZ-NRDOT) *YV- (MASS*XG-YRDOT) *NV) /DH 
A12= { (IZ-NRDOT) * (-MASS+YR) - 
& (MASS*XG-YRDOT) * ( -MASS*XG+NR) ) /DH 

A21= ( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) *YV) /DH 
A22= ( (MASS-YVDOT) * ( -MASS*XG+NR) - 
& (MASS*XG-NVDOT) * ( -MASS+YR) ) /DH 

Bll= ( (IZ-NRDOT) *YDRS- (MASS*XG-YRDOT) *NDRS) /DH 
B12= ( (IZ-NRDOT) *YDRB- (MASS*XG-YRDOT) *NDRB) /DH 
B21= ( (MASS-YVDOT) *NDRS- (MASS*XG-NVDOT) *YDRS) /DH 
B22= ( (MASS-YVDOT) *NDRB- (MASS*XG-NVDOT) *YDRB) /DH 
B1 =B11-B12 
B2 =B21-B22 



200 


WRITE 


(*,1004) 




READ 


(*,*) WNMIN,WNMAX, 




INCR=IWN 




WRITE 


( * , * ) ' ENTER OBSERVER 




READ ( 


* , * ) WNO 


205 


WRITE ( 


*,1007) 




READ ( 


*,*) A3 




DO is 


Dsat 


50 


WRITE 


(*,1006) 




READ 


(*,*) DO 




WRITE 


(*,1008) 




READ 


(*,*) CCRIT 



WN' 



Cl= (A11*A22-A21*A12) * (A21*B1-A11*B2 ) 

C2=(A11+A22) * (A21*B1-A11*B2)+B2* (A11*A22 -A21*A12 ) 

C3=- (A21*B1-A11*B2) **2 

A=C1/C2 

B=C3/C2 



DO 1 11=1, INCR 



WN =WNMIN+ (WNMAX-WNMIN) * (II-l) / (INCR-1) 
C print *,wn 

ALPHA0=WN**3 
ALPHA1=2 . 15*WN**2 
ALPHA2=1 .75*WN 
KPSI=-ALPHA1/B 
KY =-ALPHA0/B 
KR =- (ALPHA2+A) /B 
GAMAO=WNO**3 
GAMA1=2 . 15*WNO**2 
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GAMA2=1.75*WNO 

LY=A+GAMA2 

LPSI=A*LY+GAMA1 

LR=A*LPSI+GAMAO 

C234567 

C A MATRIX 

AMAT { 1 , 1 ) =0 . 0 

AMAT ( 1 , 2 ) =1 . 0 

AMAT (1, 3) =0 .0 

AMAT (1,4) =0.0 

AMAT { 1 , 5 ) =0 . 0 

AMATd, 6) =0 . 0 

AMAT(2, 1) =B*CCRIT*KPSI 

AMAT{2,2) =A+ (B*KR) 

AMAT{2, 3) =B*KY 
AMAT(2, 4) =-B*CCRIT*KPSI 
AMAT{2, 5) =-B*KR 
AMAT{2, 6) =-B*KY 
AMAT ( 3 , 1 ) =1 . 0 
AMAT (3, 2) =0.0 
AMAT (3, 3) =0 .0 
AMAT ( 3 , 4 ) =0 . 0 
AMAT ( 3 , 5 ) =0 . 0 
AMAT(3, 6)=0.0 
AMAT (4, 1) =0 .0 
AMAT (4, 2) =0.0 
AMAT (4, 3) =0 . 0 
AMAT (4, 4) =0.0 
AMAT (4, 5) =1 . 0 
AMAT (4, 6) =-LPSI 
AMAT (5, 1) =0 .0 
AMAT (5, 2) =0 .0 
AMAT ( 5 , 3 ) =0 . 0 
AMAT (5, 4) =0 . 0 
AMAT ( 5 , 5 ) =A 
AMAT(5, 6) =-LR 
AMAT(6, 1) =0 .0 
AMAT (6, 2) =0.0 
AMAT (6, 3) =0.0 
AMAT (6, 4) =1 . 0 
AMAT{6, 5) =0 . 0 
AMAT{6, 6) =-LY 
DO 11 1=1,6 
DO 12 J=l, 6 

ASAVE ( I , J) =AMAT { I , J) 

12 CONTINUE 

11 CONTINUE 

CALL RG ( 6 , 6 , AMAT, WR, WI , 1 , YYY, IVl , FVl , lERR) 
CALL DSOMEG (IEV,WR,WI, OMEGA, CHECK) 

C WRITE ( * , * ) lEV 
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5 

13 

21 

14 

22 

15 

23 

16 

24 

17 

25 

30 

C 

C 

C 



WRITE (60,*) (WR(IREAL) , IREAL=1, 6) 

OMEGAO=OMEGA 

DO 5 1=1,6 

T(I, 1) =YYY(I, lEV) 

T(I, 2) =-YYY(I, IEV+1) 

CONTINUE 

IF(IEV.EQ.l) GO TO 13 
IF(IEV.EQ.2) GO TO 14 
IF(IEV.EQ.3) GO TO 15 
IF(IEV.EQ.4) GO TO 16 
IF(IEV.EQ.5) GO TO 17 
STOP 3004 
DO 21 1=1,6 

T ( I , 3 ) =YYY (1,3) 

T ( I , 4 ) =YYY (1,4) 

T(I, 5) =YYY(I, 5) 

T ( I , 6 ) =YYY (1,6) 

CONTINUE 
GO TO 30 
DO 22 1=1,6 

T(I,3)=YYY(I, 1) 

T ( I , 4 ) =YYY (1,4) 

T(I,5)=YYY(I,5) 

T ( I , 6 ) =YYY (1,6) 

CONTINUE 
GO TO 30 
DO 23 1=1,6 

T ( I , 3 ) =YYY (1,1) 

T ( I , 4 ) =YYY (1,2) 

T(I,5)=YYY(I,5) 

T ( I , 6 ) =YYY (1,6) 

CONTINUE 
GO TO 30 
DO 24 1=1,6 

T ( I , 3 ) =YYY (1,1) 

T ( I , 4 ) =YYY (1,2) 

T ( I , 5 ) =YYY (1,3) 

T ( I , 6 ) =YYY (1,6) 

CONTINUE 
GO TO 30 
DO 25 1=1,6 

T ( I , 3 ) =YYY (1,1) 

T ( I , 4 ) =YYY (1,2) 

T ( I , 5 ) =YYY (1,3) 

T ( I , 6 ) =YYY (1,4) 

CONTINUE 

CONTINUE 

NORMALIZATION OF THE CRITICAL EIGENVECTOR 
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CALL NORMAL (T) 



C 

C INVERT TRANSFORMATION MATRIX 

C 

DO 2 1=1,6 
DO 3 J=l, 6 

TINVd, J) =0 . 0 
TSAVEd, J)=T(I, J) 

3 CONTINUE 

2 CONTINUE 

CALL DLUD (6,6, TS AVE , 6 , TLUD , I VLUD ) 

DO 4 1=1,6 

IF dVLUD(I) .EQ. 0) STOP 3003 

4 CONTINUE 

CALL DILU ( 6 , 6 , TLUD , IVLUD , SVLUD ) 

DO 8 1=1,6 
DO 9 J=l,6 

TINVd, J) =TLUD(I, J) 

9 CONTINUE 

8 CONTINUE 

C 

C CHECK Inv(T)*A*T 

C 

IMULT=1 

IF (IMULT.EQ.l) CALL MULT (TINV, ASAVE, T, A2 ) 
IF (IMULT.EQ.O) STOP 3007 
P1=A2 (1,1) 

P2=A2 (2,2) 

P=A2 (3,3) 

Q=A2 (4,4) 

R=A2 (5,5) 

S=A2 (6, 6) 

WRITE(21, *) PI, P2, P,Q,R,S 
C 

C DEFINITION OF Nij 

C 

N11=TINV(1, 1) 

N21=TINV(2, 1) 

N31=TINV(3 , 1) 

N41=TINV(4, 1) 

N51=TINV(5, 1) 

N61=TINV(6, 1) 

N12=TINV(1, 2) 

N22=TINV(2,2) 

N32=TINV(3,2) 

N42=TINV(4, 2) 

N52=TINV(5, 2) 

N62=TINV(6, 2) 

N13=TINV(1, 3) 

N23=TINV(2, 3) 
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N33=TINV(3, 3) 
N43=TINV(4, 3) 
N53=TINV{5,3) 
N63=TINV(6,3) 
N14=TINV(1, 4) 
N24=TINV{2,4) 
N34=TINV(3,4) 
N44=TINV(4, 4) 
N54=TINV(5, 4) 
N64=TINV(6, 4) 
N15=TINV(1, 5) 
N25=TINV(2, 5) 
N35=TINV(3, 5) 
N45=TINV(4, 5) 
N55=TINV(5, 5) 
N65=TINV(6, 5) 
N16=TINV(1, 6) 
N26=TINV(2, 6) 
N36=TINV(3, 6) 
N46=TINV(4, 6) 
N56=TINV(5, 6) 
N66=TINV(6, 6) 

DEFINITION OF Mij 

M11=T(1, 1) 
M21=T(2, 1) 
M31=T(3, 1) 
M41=T(4, 1) 
M51=T(5, 1) 
M61=T(6, 1) 
M12=T(1,2) 
M22=T(2,2) 
M32=T(3,2) 
M42=T(4,2) 
M52=T(5,2) 
M62=T(6,2) 
M13=T(1, 3) 
M23=T(2, 3) 
M33=T(3, 3) 
M43=T(4, 3) 
M53=T(5,3) 
M63=T(6, 3) 
M14=T(1, 4) 
M24=T(2, 4) 
M34=T(3, 4) 
M44=T(4,4) 
M54=T(5, 4) 
M64=T(6,4) 
M15=T(1, 5) 
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M25=T 


(2, 


5) 








M35=T 


(3, 


5) 








M45=T 


(4, 


5) 








M55=T 


(5, 


5) 








M65=T 


(6, 


5) 








M16=T 


(1, 


6) 








M2 6=T 


(2, 


6) 








M3 6=T 


(3, 


6) 








M46=T 


(4, 


6) 








M56=T 


(5, 


6) 








M66=T 


(6, 


6) 








WRITE 


(70 


,*) 


Nil, 


-N12, 


,N13 


WRITE 


(71 


,*) 


N14, 


-N15, 


,N16 


WRITE 


(72 


/*) 


N21, 


-N22, 


,N23 


WRITE 


(73 


,*) 


N24, 


-N25, 


,N26 


WRITE 


(74 




N31, 


-N32, 


,N33 


WRITE 


(75 


,*) 


N34, 


-N35, 


,N36 


WRITE 


(76 


,*) 


N41, 


,N42, 


,N43 


WRITE 


(77 


,*) 


N44, 


,N45, 


,N46 


WRITE 


(78 


,*) 


N51, 


-N52, 


,N53 


WRITE 


(79 


,*) 


N54, 


-N55, 


,N56 


WRITE 


(80 


,*) 


N61, 


,N62, 


,N63 


WRITE 


(81 


/ *) 


N64, 


-N65, 


,N66 



Kl=l . /8 . * ( (ALPHA2**3) +ALPHAO) / (ALPHA2) 

K2=3 . *A3- .5* (ALPHA2**2) /ALPHAO 

K3=l./{ (B**2) * (D0**2) ) * (ALPHA2+A) * ( (ALPHA0/ALPHA2 ) + (A**2 ) ) 
K=K1* (K2+K3) 

print *, wn,k,ccrit 

BL=B/ (3*D0**2) 

C2345678901234567890123456789012345678901234567890123456789012345 

6789012 

L21=A3 *M21**3-BL* (CCRIT**3*KPSI**3*M41**3+KR**3*M51**3+ 

& KY**3*M61**3+ 

& 3*CCRIT**2*KPSI**2*KR*M41**2*M51+ 

& 3*CCRIT**2*KPSI**2*KY*M41**2*M61+ 

& 3*CCRIT*KPSI*KR**2*M51**2*M41+3*CCRIT*KPSI*KY**2*M61**2*M41+ 
Sc 3*KR*KY**2*M61**2*M51+3*KR**2*KY*M51**2*M61+ 

Sc 6*CCRIT*KPSI*KR*KY*M41*M51*M61) 

L22=A3*M22**3-BL* (CCRIT**3*KPSI**3*M42**3+KR**3*M52**3+ 

Sc KY**3*M62**3 + 

Sc 3*CCRIT**2*KPSI**2*KR*M42**2*M52 + 

Sc 3*CCRIT**2*KPSI**2*KY*M42**2*M62 + 

& 3 *CCRIT*KPSI *KR* *2 *M52 * *2 *M42+3 *CCRIT*KPSI*KY* * 2 *M62 * *2 *M42 + 
Sc 3*KR*KY**2*M62**2*M52+3*KR**2*KY*M52**2*M62+ 

Sc 6*CCRIT*KPSI*KR*KY*M42*M52*M62) 

L23=3*A3*M21**2*M22-BL* (3*CCRIT**3*KPSI**3*M41**2*M42+ 

Sc 3*KR**3*M51**2*M52 + 
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& 3*KY**3*M61**2*M62+ 

& 3*CCRIT**2*KPSI**2*KR* (M41 * *2 *M52+2 *M41 *M42 *M51 ) + 

Sc 3*CCRIT**2*KPSI**2*KY* (M41**2*M62+2*M41*M42*M61 ) + 

& 3*CCRIT*KPSI*KR**2* (M51**2*M42+2*M51*M52*M41) + 

& 3*CCRIT*KPSI*KY**2* (M61**2*M42+2*M61*M62*M41 ) + 

Sc 3*KR*KY**2* (M61**2*M52+2*M61*M62*M51) + 

Sc 3*KR**2*KY* (M51**2*M62+2*M51*M52*M61) + 

Sc 6*CCRIT*KPSI**KR*KY* (M41*M51*M62+ (M41*M52+M42*M51 ) *M61) ) 
L24=3*A3*M21*M22**2-BL* (3*CCRIT**3*KPSI**3*M41*M42**2+ 

Sc 3*KR**3*M51*M52**2 + 

Sc 3*KY**3*M61*M62**2 + 

& 3*CCRIT**2*KPSI**2*KR* (M42 **2*M51+2 *M41*M42 *M52 ) + 

Sc 3*CCRIT**2*KPSI**2*KY* (M42**2*M61+2*M41*M42*M62 ) + 

Sc 3*CCRIT*KPSI*KR**2* (M52**2*M41+2*M51*M52*M42 ) + 

Sc 3*CCRIT*KPSI*KY**2* (M62**2*M41+2*M61*M62*M42 ) + 

& 3*KR*KY**2* (M62**2*M51+2*M61*M62*M52) + 

Sc 3*KR**2*KY* (M52**2*M61+2*M51*M52*M62) + 

Sc 6*CCRIT*KPSI*KR*KY* (M42*M52*M61+ (M41*M52+M42*M51) *M62) ) 

L31=(-l./6. ) *M11**3 
L32=(-l./6. ) *M12**3 
L33=(-l./2 . ) *M11**2*M12 
L34=(-l./2. ) *M11*M12**2 
R11=N12*L21+N13*L31 
R12=N12*L23+N13*L33 
R13=N12*L24+N13*L34 
R14=N12*L22+N13*L32 
R21=N22*L21+N23*L31 
R22=N22*L23+N23*L33 
R23=N22*L24+N23*L34 
R24=N22*L22+N23*L32 



C 



K= (3*R11+R13+R22+3*R24) /8 
WRITE (11,2001) WN,K,CCRIT 
1 CONTINUE 
STOP 

1004 FORMAT (' ENTER MIN, MAX, AND INCREMENTS OF 



1006 FORMAT 

1007 FORMAT 

1008 FORMAT 
2001 FORMAT 

END 



WN(stepstodivide) ' ) 
( ' ENTER DELTASAT . ' ) 

( ' ENTER A3 ' ) 

(' ENTER CCRIT') 
(3E15.5) 



SUBROUTINE DSOMEG ( I JK , WR , WI , OMEGA , CHECK ) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DIMENSION WR ( 6 ) , WI ( 6 ) 

CHECK=-1 . OD+25 
DO 1 1=1,6 

IF (WR(I) .LT. CHECK) GO TO 1 
CHECK=WR ( I ) 
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IJ=I 

1 CONTINUE 

OMEGA=DABS (WI (IJ) ) 

IF (WI (IJ) .GT. 0 .DO) IJK=IJ 



IF (WI (IJ) .LT.O.DO) IJK=IJ-1 

RETURN 

END 

C====================================================== 

SUBROUTINE NORMAL (T) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION T ( 6 , 6 ) , TNOR (6,6) 

CFAC=T(1, 1) **2+T(l, 2) **2 

IF (DABS (CFAC) .LE. (1 .D-10) ) STOP 4001 

TNOR(l,l)=l.D0 

TNOR (2,1)=(T(1,1)*T(2,1)+T(2,2)*T(1,2)) /CFAC 
TNOR (3, 1) = (T(l, 1) *T(3, 1) +T(3, 2) *T(1, 2) ) /CFAC 
TNOR (4, 1) = (T(l, 1) *T(4, 1) +T(4, 2) *T(1, 2) ) /CFAC 
TNOR (5,1)=(T(1,1)*T(5,1)+T(5,2)*T(1,2)) /CFAC 
TNOR (6, 1) = (T(1,1)*T(6,1)+T(6,2)*T(1,2)) /CFAC 
TNOR(l,2) =0.D0 

TNOR(2,2) =(T(2,2)*T(1,1) -T(2,l) *T(1,2) ) /CFAC 
TNOR(3,2)=(T(3,2) *T(1,1) -T(3,l) *T(1,2) ) /CFAC 
TNOR(4,2) = (T(4,2)*T(1,1)-T(4,1)*T(1,2)) /CFAC 
TNOR(5,2)=(T(5,2) *T(1,1) -T(5,l) *T(1,2) ) /CFAC 
TNOR (6,2)=(T(6,2)*T(1,1)-T(6,1)*T(1,2)) /CFAC 
DO 1 1=1,6 
DO 2 J=l,2 

T(I, J) =TNOR(I, J) 

2 CONTINUE 

1 CONTINUE 
RETURN 
END 

C=== ========== ========================================= 

SUBROUTINE MULT (TINV, A, T, A2 ) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION TINV (6,6),A(6,6),T(6,6),A1(6,6),A2(6,6) 
DO 1 1=1,6 
DO 2 J=l, 6 
A1 (I, J) =0 .DO 
A2 (I, J) =0 .DO 

2 CONTINUE 

1 CONTINUE 

DO 3 1=1,6 
DO 4 J=l,6 
DO 5 K=l,6 

A1 (I, J) =A(I,K) *T(K, J) +A1 (I, J) 
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LD ^ ro 



CONTINUE 
CONTINUE 
CONTINUE 
DO 6 1=1,6 
DO 7 J=l,6 
DO 8 K=l,6 

A2 (I, J) =TINV(I,K) *A1 (K, J) +A2 (I, J) 
8 CONTINUE 



7 CONTINUE 
6 CONTINUE 

DO 11 1=1,6 



c 


WRITE (*,101) 


(A(I, J) , J=l, 6) 


11 


CONTINUE 






DO 12 1=1,6 




c 


WRITE (*,101) 


(T(I, J) , J=l, 6) 


12 


CONTINUE 






DO 10 1=1,6 




C 


WRITE (*,101) 


(A2 (I, J) , J=l, 6) 


10 


CONTINUE 




C 


WRITE (*,101) 


A2 ( 1 , 1 ) 




RETURN 




101 


FORMAT (4E15 


.5) 




END 
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APPENDIX B 

CRITICAL VALUE OF C FOR [X, X] BASIS 



PROGRAM NCCRIT 
HOPF BIFURCATIONS 
NOMOTO'S FIRST ORDER MODEL 

CALCULATIONS FOR CCRITICAL 

C234567 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION Kl , K2 , K, LPSI , LY, LR, IZ , L, 

& MASS,NV,NR,NVDOT,NRDOT,NDRS,NDRB, KPSI, KR, KY, K3 

C 

DIMENSION AMAT (6,6), FVl ( 6 ) , IVl ( 6 ) , YYY (6,6) 

DIMENSION WR(6) ,WI(6) ,TSAVE(6,6) ,TLUD(6,6) ,IVLUD(6) ,SVLUD(6) 
DIMENSION ASAVE ( 6 , 6 ) , Al ( 6 , 6 ) , A2 ( 6 , 6 ) 

OPEN (11,FILE='CVWN1.RES' , STATUS =' NEW ' ) 

OPEN (12,FILE='CVWN2 .RES' ,STATUS='NEW' ) 

OPEN (13,FILE='CVWN3.RES' , STATUS =' NEW ' ) 

WRITE (*,*) 'ENTER MIN, MAX, AND INCREMENTS IN CCRIT' 

READ (*,♦) CMIN,CMAX,IC 

WRITE (*,*) 'ENTER MIN, MAX, AND INCREMENTS IN WN' 

READ (*,*) WNMIN,WNMAX, INCR 
WRITE ( * , ♦ ) ' ENTER WNO ' 

READ ( * , ♦ ) WNO 

WEIGHT=435.0 

IZ =45.0 

L =7.3 

RHO =1.94 

G =32.2 

XG =0.0104 

MASS =WEIGHT/G 

MASS =MASS/ (0.5*RHO*L**3) 

IZ =IZ/ (0 .5*RHO*L**5) 

XG =XG/L 
YRDOT =-0.00000 
YVDOT =-0.03430 
YR =+0.00000 
YV =-0.10700 
YDRS =+0.01241 
YDRB =+0.01241 
NRDOT =-0.00047 
NVDOT =-0.00000 
NR =-0.00390 
NV =-0.00000 
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NDRS =-0.337*YDRS 
NDRB =+0.283*YDRS 
DH = (IZ-NRDOT) * (MASS-YVDOT) - 
& (MASS*XG-YRDOT) * (MASS*XG-NVDOT) 

All= ( (IZ-NRDOT) *YV- (MASS*XG-YRDOT) *NV) /DH 



A12= ( (IZ-NRDOT) * (-MASS+YR) - 
& (MASS*XG-YRDOT) * ( -MASS*XG+NR) ) /DH 

A21= ( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) *YV) /DH 
A22= ( (MASS-YVDOT) * ( -MASS*XG+NR) - 
& (MASS*XG-NVDOT) * ( -MASS+YR) ) /DH 

Bll= ( (IZ-NRDOT) *YDRS- (MASS*XG-YRDOT) *NDRS) /DH 
B12= ( (IZ-NRDOT) *YDRB- (MASS*XG-YRDOT) *NDRB) /DH 
B21= ( (MASS-YVDOT) *NDRS- (MASS*XG-NVDOT) *YDRS) /DH 
B22= ( (MASS-YVDOT) *NDRB- (MASS*XG-NVDOT) *YDRB) /DH 
B1 =B11-B12 
B2 =B21-B22 

Cl= (A11*A22-A21*A12) * ( A21*B1-A11*B2 ) 

C2= (A11+A22) * (A21*B1-A11*B2) +B2* (A11*A22-A21*A12 ) 
C3=- (A21*B1-A11*B2) **2 
A=C1/C2 
B=C3/C2 
C 

EPS=1 .D-5 
ILMAX=1500 
C 

DO 1 11=1, INCR 
C 

WN =WNMIN+ (WNMAX-WNMIN) * (II-l) / (INCR-1) 

C print *,wn 

ALPHA0=WN**3 
ALPHA1=2 .15*WN**2 
ALPHA2=1.75*WN 
KPSI=-ALPHA1/B 
KY =-ALPHA0/B 
KR =- (ALPHA2+A) /B 
GAMAO=WNO**3 
GAMA1=2 . 15*WNO**2 
GAMA2=1.75*WNO 
LY=A+GAMA2 
LPSI=A*LY+GAMA1 
LR=A*LPSI+GAMAO 
C234567 

DO 2 J=1,IC 

CCRIT=CMIN+ ( J-1) * (CMAX-CMIN) / (IC-1) 

C A MATRIX 

AMATd, 1) =0 . 0 
AMATd, 2) =1 . 0 
AMATd, 3) =0 .0 
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AMATd, 4) =0.0 

AMATd, 5) =0.0 

AMATd, 6) =0 . 0 

AMAT(2,1)=0.0 

AMAT ( 2 , 2 ) =A 

AMAT(2,3)=0.0 

AMAT(2 , 4) =B*CCRIT*KPSI 

AMAT (2, 5) =B*KR 

AMAT(2, 6) =B*KY 

AMAT { 3 , 1 ) =1 . 0 

AMAT (3, 2) =0.0 

AMAT(3,3)=0.0 

AMAT (3, 4) =0.0 

AMAT(3,5)=0.0 

AMAT (3 , 6) =0 . 0 

AMAT (4, 1) =0 .0 

AMAT (4, 2) =0.0 

AMAT(4, 3) =LPSI 

AMAT (4, 4) =0.0 

AMAT (4, 5) =1.0 

AMAT(4, 6) =-LPSI 

AMAT (5, 1) =0.0 

AMAT (5, 2) =0 .0 

AMAT (5,3) =LR 

AMAT (5,4) =B*CCRIT*KPSI 

AMAT(5, 5) =B*KR 

AMAT (5,6) =-LR+B*KY 

AMAT (6, 1) =0.0 

AMAT (6, 2) =0.0 

AMAT (6,3) =LY 

AMAT (6, 4) =1 . 0 

AMAT (6, 5) =0 .0 

AMAT(6, 6) =-LY 

C 

CALL RG ( 6 , 6 , AMAT, WR, WI , 0 , ZZZ, IVl , FVl , lERR) 
CALL DSTABL(DEOS,WR,WI,FREQ) 

U=CCRIT 

IF (J.GT.l) GO TO 10 

DEOSOO=DEOS 

UOO=U 

LL=0 

GO TO 2 

10 DEOSNN=DEOS 
UNN=U 

PR=DEOSNN*DEOSOO 

IF (PR.GT.O.DO) GO TO 3 

LL=LL+1 

IF (LL.GT.3) STOP 1000 

IL=0 

UO=UOO 
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UN=UNN 
DEOSO=DEOSOO 
DE0SN=DE0S1SIN 
6 UL=UO 

UR=UN 

DEOSL=DEOSO 
DEOSR=DEOSN 
U= (UL+UR) /2 .DO 
CCRIT=U 
AMATd, 1) =0 . 0 
AMAT ( 1 , 2 ) =1 . 0 
AMAT ( 1 , 3 ) =0 . 0 
AMAT ( 1 , 4 ) =0 . 0 
AMATd, 5) =0 .0 
AMATd, 6) =0 .0 
AMAT (2, 1) =0.0 
AMAT ( 2 , 2 ) =A 
AMAT (2, 3) =0.0 
AMAT (2,4) =B*CCRIT*KPSI 
AMAT(2, 5) =B*KR 
AMAT (2, 6) =B*KY 
AMAT ( 3 , 1 ) =1 . 0 
AMAT (3, 2) =0.0 
AMAT (3, 3) =0 .0 
AMAT (3, 4) =0.0 
AMAT (3, 5) =0 .0 
AMAT (3, 6) =0 .0 
AMAT(4,1)=0.0 
AMAT ( 4 , 2 ) =0 . 0 
AMAT(4, 3) =LPSI 
AMAT (4, 4) =0 .0 
AMAT ( 4 , 5 ) =1 . 0 
AMAT(4, 6) =-LPSI 
AMAT ( 5 , 1 ) =0 . 0 
AMAT ( 5 , 2 ) =0 . 0 
AMAT (5,3) =LR 
AMAT (5,4) =B*CCRIT*KPSI 
AMAT(5, 5) =B*KR 
AMAT(5, 6) =-LR+B*KY 
AMAT ( 6 , 1 ) =0 . 0 
AMAT (6, 2) =0 . 0 
AMAT (6,3) =LY 
AMAT (6, 4) =1 . 0 
AMAT ( 6 , 5 ) =0 . 0 
AMAT (6, 6) =-LY 
C 

CALL RG(6, 6, AMAT,WR,WI, 0, ZZZ, IVl , FVl , lERR) 
CALL DSTABL ( DECS, WR,WI, FREQ) 

U=CCRIT 

DEOSM=DEOS 
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UM=U 

PRL=DEOSL*DEOSM 

PRR=DEOSR*DEOSM 

IF (PRL.GT.O.DO) GO TO 5 

UO=UL 

UN=UM 

DEOSO=DEOSL 

DEOSN=DEOSM 

IL=IL+1 

IF (IL.GT.ILMAX) STOP 3100 
DIF=DABS (UL-UM) 

IF (DIF.GT.EPS) GO TO 6 

U=UM 

GO TO 4 

5 IF (PRR.GT.O.DO) STOP 3200 

UO=UM 
UN=UR 

DEOSO=DEOSM 

DEOSN=DEOSR 

IL=IL+1 

IF (IL.GT.ILMAX) STOP 3100 
DIF=DABS (UM-UR) 

IF (DIF.GT.EPS) GO TO 6 
TJ=UM 

4 LLL=10+LL 

CCRIT=U 

WRITE (LLL,*) CCRIT,WN 
3 UOO=UNN 

DEOSOO=DEOSNN 
2 CONTINUE 

1 CONTINUE 

C 

2001 FORMAT (215) 

END 



SUBROUTINE DSTABL (DEOS , WR, WI , OMEGA) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DIMENSION WR ( 6 ) , WI ( 6 ) 

DEOS=-1.0D+20 
DO 1 1=1,6 

IF (WR(I) .LT.DEOS) GO TO 1 
DEOS=WR(I) 

IJ=I 

1 CONTINUE 

OMEGA=WI (IJ) 

OMEGA=DABS (OMEGA) 

RETURN 

END 

C== ====================================== 
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APPENDIX C 

HOPF BIFURCATION PROGRAM FOR [X,X] BASIS 



PROGRAM KKPSI 

HOPF BIFURCATIONS 
NOMOTO'S FIRST ORDER MODEL 

CALCULATIONS FOR K AND CCRITICAL IF KPSI CHANGES W/ C 



234567 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 



DOUBLE PRECISION K1 , K2 , K, LPSI , LY, LR, IZ , L, 

& MASS , NV, NR, NVDOT, NRDOT, NDRS , NDRB, KPSI , KR, KY, K3 , 

& L21,L22,L23,L24,L31,L32,L33,L34,L51,L52,L53,L54, 

& M11,M12,M13,M14,M15,M16,M21,M22,M23,M24,M25,M26, 

& M31 , M32 , M33 , M34 , M3 5 , M3 6 , M41 , M42 , M43 , M44 , M45 , M46 , 

& M51,M52,M53,M54,M55,M56,M61,M62,M63,M64,M65,M66, 

& N11,N12,N13,N14,N15,N16,N21,N22,N23,N24,N25,N26, 

& N3 1 , N3 2 , N3 3 , N3 4 , N3 5 , N3 6 , N4 1 , N4 2 , N4 3 , N4 4 , N4 5 , N4 6 , 

& N51 , N52 , N53 , N54 , N55 , N56 , N61 , N62 , N63 , N64 , N65 , N66 



DIMENSION AMAT (6,6),T{6,6), TINV (6,6), FVl ( 6 ) , IVl ( 6 ) , YYY (6,6) 
DIMENSION WR ( 6 ) , WI ( 6 ) , TSAVE (6,6), TLUD (6,6), IVLUD ( 6 ) , SVLUD ( 6 ) 
DIMENSION ASAVE ( 6 , 6 ) , Al ( 6 , 6 ) , A2 ( 6 , 6 ) 

OPEN dO,FILE='CVWNl.RES' ,STATUS='OLD' ) 

OPEN (11,FILE='AKPSI.MAT' , STATUS= 'NEW' ) 

OPEN (12,FILE='IREAL.MAT' , STATUS =' NEW' ) 

OPEN (13,FILE='RVALS.MAT' ,STATUS='NEW' ) 

OPEN (14,FILE='KKKKK.MAT' , STATUS= ' NEW' ) 



WEIGHT=435. 0 


IZ 


=45.0 


L 


=7.3 


RHO 


= 1.94 


G 


=32.2 


:g 


=0.0104 


i-lASS 


=WEIGHT/G 


MASS 


=MASS/ (0.5*RHO*L** 


IZ 


=IZ/ (0 .5*RHO*L**5) 


XG 


=XG/L 


YRDOT 


=-0.00000 


YVDOT 


=-0.03430 


YR 


=+0.00000 


YV 


=-0.10700 


YDRS 


=+0.01241 


YDRB 


=+0.01241 
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NRDOT =-0.00047 
NVDOT =-0.00000 
NR =-0.00390 

NV =-0.00000 

NDRS =-0.337*YDRS 
NDRB =+0.283*YDRS 
DH = (IZ-NRDOT) * (MASS-YVDOT) - 
& (MASS*XG-YRDOT) * ( MAS S*XG- NVDOT) 

All={ (IZ-NRDOT) *YV- (MASS*XG-YRDOT) *NV) /DH 
A12= { (IZ-NRDOT) * (-MASS+YR) - 
& (MASS*XG-YRDOT) * ( -MASS*XG+NR) ) /DH 

A21= ( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) *YV) /DH 
A22= ( (MASS-YVDOT) * ( -MASS*XG+NR) - 
& ( MAS S*XG -NVDOT) * ( -MASS+YR) ) /DH 

Bll= ( (IZ-NRDOT) *YDRS- (MASS*XG-YRDOT) *NDRS) /DH 
B12= ( (IZ-NRDOT) *YDRB- (MASS*XG-YRDOT) *NDRB) /DH 
B21= ( (MASS-YVDOT) *NDRS- ( MAS S*XG -NVDOT) *YDRS) /DH 
B22= ( (MASS-YVDOT) *NDRB- ( MAS S*XG- NVDOT) *YDRB) /DH 
B1 =B11-B12 
B2 =B21-B22 
C 

200 WRITE (*,1004) 

READ ( * , * ) IWN 
INCR=IWN 

WRITE (*,*) 'ENTER OBSERVER WN' 

READ ( * , * ) WNO 
205 WRITE(*,1007) 

READ (*,*) A3 
C DO is Dsat 

50 WRITE (*,1006) 

READ ( * , * ) DO 

C1=(A11*A22-A21*A12) * (A21*Bl-All*B2 ) 

C2= (A11+A22 ) * (A21*B1-A11*B2 ) +B2* (A11*A22-A21*A12 ) 
C3=- (A21*B1-A11*B2) **2 
A=C1/C2 
B=C3/C2 
A3=0.0 



DO 1 11=1, INCR 



READ (10,*) CCRIT,WN 
C print *,wn 

ALPHA0=WN**3 
ALPHA1=2 .15*WN**2 
ALPHA2=1 .75*WN 
KPSI=-ALPHA1/B 
KY =-ALPHA0/B 
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KR =- (ALPHA2+A) /B 

GAMAO=WNO**3 

GAMA1=2 . 15*WNO**2 

GAMA2=1.75*WNO 

LY=A+GAMA2 

LPSI=A*LY+GAMA1 

LR=A*LPSI+GAMAO 

C234567 

C A MATRIX 

AMATd, 1) =0 . 0 
AMAT ( 1 , 2 ) =1 . 0 
AMAT { 1 , 3 ) =0 . 0 
AMAT ( 1 , 4 ) =0 . 0 
AMAT ( 1 , 5 ) =0 . 0 
AMATd, 6) =0 .0 
AMAT (2,1) =0.0 
AMAT { 2 , 2 ) =A 
AMAT (2, 3) =0.0 
AMAT(2, 4) =B*CCRIT*KPSI 
AMAT{2, 5) =B*KR 
AMAT (2 )=B*KY 
AMAT (2 =1.0 

AMAT (3,-; =0.0 
AMAT { 3 , 3 ) =0 . 0 
AMAT (3, 4) =0.0 
AMAT (3, 5) =0 .0 
AMAT ( 3 , 6 ) =0 . 0 
AMAT(4,1)=0.0 
AMAT (4, 2) =0.0 
AMAT(4, 3) =LPSI 
AMAT ( 4 , 4 ) =0 . 0 
AMAT(4,5)=1.0 
AMAT (4, 6) =-LPSI 
AMAT(5,1)=0.0 
AMAT ( 5 , 2 ) =0 . 0 
AMAT (5,3) =LR 
AMAT (5,4) =B*CCRIT*KPSI 
AMAT (5, 5) =B*KR 
AMAT(5, 6) =-LR+B*KY 
AMAT (6, 1) =0 .0 
AMAT ( 6 , 2 ) =0 . 0 
AMAT (6,3) =LY 
AMAT ( 6 , 4 ) =1 . 0 
AMAT (6, 5) =0.0 
AMAT ( 6 , 6 ) =-LY 
DO 11 1=1,6 
DO 12 J=l, 6 

ASAVEd, J) =AMAT(I, J) 

12 CONTINUE 

II CONTINUE 
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c 

5 

13 

21 

14 

22 

15 

23 

16 

24 

17 

25 

30 



CALL RG(6, 6, AMAT,WR,WI, 1, YYY, IVl, FVl, lERR) 
CALL DSOMEG ( lEV , WR , WI , OMEGA , CHECK ) 

WRITE (*,*) lEV 

WRITE (12,*) (WR(IREAL) , IREAL=1, 6) 

OMEGAO=OMEGA 

DO 5 1=1,6 

T(I, 1) =YYY(I, lEV) 

T ( I , 2 ) =-YYY ( I , IEV+1 ) 

CONTINUE 

IF(IEV.EQ.l) GO TO 13 
IF(IEV.EQ.2) GO TO 14 
IF(IEV.EQ.3) GO TO 15 
IF(IEV.EQ.4) GO TO 16 
IF(IEV.EQ.5) GO TO 17 
STOP 3004 
DO 21 1=1,6 

T(I,3)=YYY(I,3) 

T ( I , 4 ) =YYY (1,4) 

T ( 1 , 5 ) =YYY (1,5) 

T ( I , 6 ) =YYY (1,6) 

CONTINUE 
GO TO 30 
DO 22 1=1,6 

T(I,3)=YYY(I,1) 

T ( I , 4 ) =YYY (1,4) 

T ( I , 5 ) =YYY (1,5) 

T ( I , 6 ) =YYY (1,6) 

CONTINUE 
GO TO 30 
DO 23 1=1,6 

T ( I , 3 ) =YYY (1,1) 

T ( I , 4 ) =YYY (1,2) 

T ( 1 , 5 ) =YYY (1,5) 

T ( I , 6 ) =YYY (1,6) 

CONTINUE 
GO TO 30 
DO 24 1=1,6 

T ( 1 , 3 ) =YYY (1,1) 

T ( I , 4 ) =YYY (1,2) 

T ( I , 5 ) =YYY (1,3) 

T ( I , 6 ) =YYY (1,6) 

CONTINUE 
GO TO 30 
DO 25 1=1,6 

T ( I , 3 ) =YYY (1,1) 

T ( I , 4 ) =YYY (1,2) 

T ( I , 5 ) =YYY (1,3) 

T ( I , 6 ) =YYY (1,4) 

CONTINUE 

CONTINUE 
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c 

C NORMALIZATION OF THE CRITICAL EIGENVECTOR 

C 

CALL NORMAL (T) 

C 

C INVERT TRANSFORMATION MATRIX 

C 

DO 2 1=1,6 
DO 3 J=l, 6 

TINVd, J)=0.0 
TSAVEd, J) =T(I, J) 

3 CONTINUE 

2 CONTINUE 

CALL DLUD (6,6, TSAVE , 6 , TLUD , IVLUD) 

DO 4 1=1, 6 

IF (IVLUDd) .EQ.O) STOP 3003 

4 CONTINUE 

CALL DILU ( 6 , 6 , TLUD , IVLUD , SVLUD ) 

DO 8 1=1, 6 
DO 9 J=l, 6 

TINV (I , J) =TLUD (I , J ) 

9 CONTINUE 

8 CONTINUE 

C 

C CHECK Inv(T)*A*T 

C 

IMULT=1 

IF (IMULT.EQ.l) CALL MULT (TINV, AS AVE, T, A2 ) 
IF (IMULT.EQ.O) STOP 
P=A2 (3,3) 

Q=A2(4,4) 

R=A2 (5,5) 

S=A2 (6,6) 

C WRITE(21,*) P,Q,R,S 

C DEFINITION OF Nij 

C 

N11=TINV(1, 1) 

N21=TINV(2, 1) 

N31=TINV(3, 1) 

N41=TINV(4, 1) 

N51=TINV(5, 1) 

N61=TINV(6, 1) 

N12=TINV(1, 2) 

N22=TINV(2, 2) 

N32=TINV(3, 2) 

N42=TINV(4, 2) 

N52=TINV(5, 2) 

N62=TINV(6, 2) 

N13=TINV(1, 3) 
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N23=TINV(2, 3) 
N33=TINV(3 , 3 ) 
N43=TINV(4, 3) 
N53=TINV(5, 3) 
N63=TINV(6, 3) 
N14=TINV(1, 4) 
N24=TINV(2, 4) 
N34=TINV(3 , 4) 
N44=TINV(4, 4) 
N54=TINV(5, 4) 
N64=TINV(6, 4) 
N15=TINV(1, 5) 
N25=TINV(2, 5) 
N35=TINV(3, 5) 
N45=TINV(4, 5) 
N55=TINV(5, 5) 
N65=TINV(6, 5) 
N16=TINV(1, 6) 
N26=TINV(2, 6) 
N36=TINV(3 , 6) 
N46=TINV(4, 6) 
N56=TINV(5, 6) 
N66=TINV(6, 6) 

DEFINITION OF Mij 

M11=T(1, 1) 
M21=T{2, 1) 
M31=T(3, 1) 
M41=T{4, 1) 
M51=T(5, 1) 
M61=T{6, 1) 
M12=T(1,2) 
M22=T(2,2) 
M32=T(3,2) 
M42=T(4, 2) 
M52=T(5, 2) 
M62=T(6,2) 
M13=T(1, 3) 
M23=T(2, 3) 
M33=T(3, 3) 
M43=T(4,3) 
M53=T(5, 3) 
M63=T(6, 3) 
M14=T(1, 4) 

M24=T(2 , 4) 

M34=T(3 , 4) 
M44=T(4, 4) 
M54=T(5, 4) 
M64=T(6, 4) 
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M15=T(1, 5) 

M25=T(2, 5) 

M35=T(3, 5) 

M45=T(4, 5) 

M55=T(5, 5) 

M65=T(6, 5) 

M16=T(1, 6) 

M26=T(2, 6) 

M36=T(3, 6) 

M46=T(4, 6) 

M56=T{-j, 6) 

M66=T ^ 6, 6) 

Kl=l. /8. * ( (ALPHA2**3) +ALPHA0) / (ALPHA2) 

K2=3 .*A3-.5* (ALPHA2**2) /ALPHAO 

K3=l . / ( (B**2) * (D0**2) ) * (ALPHA2+A) * ( ( ALPHAO /ALPHA2 ) + (A**2) ) 
K=K1* (K2+K3) 

print *, wn,k,ccrit 

BL=B/ (3*D0**2) 

C 

C2345678901234567890123456789012345678901234567890123456789012345 
L21=A3*M21**3-BL* (CCRIT**3*KPSI**3*M41**3+KR**3 *M51**3+ 

& KY**3*M61**3+ 

& 3*CCRIT**2*KPSI**2*KR*M41**2*M51+ 

& 3*CCRIT**2*KPSI**2*KY*M41**2*M61+ 

Sc 3*CCRIT*KPSI*KR**2*M51**2*M41+3*CCRIT*KPSI*KY**2*M61**2*M41+ 
Sc 3*KR*KY**2*M61**2*M51+3*KR**2*KY*M51**2*M61+ 

Sc 6*CCRIT*KPSI*KR*KY*M41*M51*M61) 

L22=A3*M22**3-BL* (CCRIT**3*KPSI**3*M42**3+KR**3*M52**3+ 

Sc KY**3*M62**3 + 

Sc 3*CCRIT**2*KPSI**2*KR*M42**2*M52 + 

Sc 3*CCRIT**2*KPSI**2*KY*M42**2*M62 + 

Sc 3 *CCRIT*KPSI *KR* *2 *M52 * *2 *M42+3 *CCRIT*KPSI*KY* *2 *M62 * *2 *M42 + 
Sc 3*KR*KY**2*M62**2*M52+3*KR**2*KY*M52**2*M62+ 

Sc 6*CCRIT*KPSI*KR*KY*M42*M52*M62) 

L23=3*A3*M21**2*M22-BL* (3*CCRIT**3*KPSI**3*M41**2*M42+ 

Sc 3*KR**3*M51**2*M52 + 

Sc 3*KY**3*M61**2*M62 + 

Sc 3*CCRIT**2*KPSI**2*KR* (M41**2*M52+2*M41*M42*M51) + 

Sc 3*CCRIT**2*KPSI**2*KY* (M41**2*M62+2*M41*M42*M61) + 

Sc 3*CCRIT*KPSI*KR**2* (M51* *2 *M42+2 *M51*M52 *M41 ) + 

Sc 3*CCRIT*KPSI*KY**2* (M61**2*M42+2*M61*M62*M41) + 

Sc 3*KR*KY**2* (M61**2*M52+2*M61*M62*M51) + 

Sc 3*KR**2*KY* (M51**2*M62+2*M51*M52*M61) + 

Sc 6*CCRIT*KPSI**KR*KY* (M41*M51*M62+ (M41*M52+M42 *M51 ) *M61) ) 

L24=3*A3*M21*M22**2-BL* (3 *CCRIT**3*KPSI**3*M41*M42**2+ 

Sc 3*KR**3*M51*M52**2 + 

Sc 3*KY**3*M61*M62**2 + 
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& 3*CCRIT**2*KPSI**2*KR* (M42 * *2 *M51+2 *M41 *M42 *M52 ) + 

& 3*CCRIT**2*KPSI**2*KY* (M42 **2*M61+2*M41*M42*M62 ) + 

& 3*CCRIT*KPSI*KR**2* {M52**2*M41+2*M51*M52*M42 ) + 

& 3*CCRIT*KPSI*KY**2* (M62**2*M41+2*M61*M62*M42 ) + 

& 3*KR*KY**2* (M62**2*M51+2*M61*M62*M52)+ 

& 3*KR**2*KY* (M52**2*M61+2*M51*M52*M62) + 

& 6*CCRIT*KPSI*KR*KY* {M42 *M52 *M61+ (M41*M52+M42 *M51 ) *M62) ) 

L31= (-1/6) *M11**3 
L32={-l/6) *M12**3 
L33= (-1/2) *M11**2*M12 
L34={-l/2) *M11*M12**2 

C2 345678901234567890123456789012345678901234567890123456789012345 
L51=-BL* (CCRIT**3*KPSI**3*M41**3+KR**3*M51**3+KY**3*M61**3+ 
Sc 3*CCRIT**2*KPSI**2*KR*M41**2*M51 + 

Sc 3*CCRIT**2*KPSI**2*KY*M41**2*M61 + 

Sc 3*CCRIT*KPSI*KR**2*M51**2*M41+3*CCRIT*KPSI*KY**2*M61**2*M41+ 
Sc 3*KR*KY**2*M61**2*M51+3*KR**2*KY*M51**2*M61 + 

& 6*CCRIT*KPSI*KR*KY*M41*M51*M61 ) 

L52=-BL* (CCRIT**3*KPSI**3*M42**3+KR**3*M52**3+KY**3*M62**3+ 
Sc 3*CCRIT**2*KPSI**2*KR*M42**2*M52 + 

Sc 3*CCRIT**2*KPSI**2*KY*M42**2*M62 + 

Sc 3*CCRIT*KPSI*KR**2*M52**2*M42+3*CCRIT*KPSI*KY**2*M62**2*M42 + 
& 3*KR*KY**2*M62**2*M52+3*KR**2*KY*M52**2*M62+ 

Sc 6*CCRIT*KPSI*KR*KY*M42*M52*M62) 

L53=-BL* (3*CCRIT**3*KPSI**3*M41**2*M42+3*KR**3*M51**2*M52+ 

Sc 3*KY**3*M61**2*M62 + 

Sc 3*CCRIT**2*KPSI**2*KR* (M41**2*M52+2*M41*M42*M51 ) + 

Sc 3*CCRIT**2*KPSI**2*KY* {M41**2*M62+2*M41*M42*M61) + 

Sc 3*CCRIT*KPSI*KR**2* {M51**2*M42+2*M51*M52*M41) + 

Sc 3*CCRIT*KPSI*KY**2* (M61**2 *M42 + 2 *M61*M62 *M41 ) + 

Sc 3*KR*KY**2* (M61**2*M52+2*M61*M62*M51) + 

Sc 3*KR**2*KY* {M51**2*M62+2*M51*M52*M61) + 

Sc 6*CCRIT*KPSI*KR*KY* (M41*M51*M62+ (M41*M52+M42 *M51 ) *M61 ) ) 
L54=-BL* (3*CCRIT**3*KPSI**3*M41*M42**2+3*KR**3*M51*M52**2+ 

Sc 3*KY**3*M61*M62**2 + 

Sc 3*CCRIT**2*KPSI**2*KR* (M42**2*M51+2*M41*M42*M52 ) + 

& 3*CCRIT**2*KPSI**2*KY* (M42**2*M61+2*M41*M42*M62 ) + 

Sc 3*CCRIT*KPSI*KR**2* (M52**2*M41+2*M51*M52*M42 ) + 

& 3*CCRIT*KPSI*KY**2* (M62**2*M41+2*M61*M62*M42 ) + 

Sc 3*KR*KY**2* (M62**2*M51+2+M61*M62*M52) + 

Sc 3*KR**2*KY* (M52**2*M61+2*M51*M52*M62) + 

Sc 6*CCRIT*KPSI*KR*KY* (M42*M52*M61+ (M41*M52+M42 *M51 ) *M62 ) ) 

R11=N12*L21+N13*L31+N15*L51 
R12=N12*L23+N13*L33+N15*L53 
R13=N12*L24+N13*L34+N15*L54 
R14=N12*L22+N13*L32+N15*L52 
R21=N22*L21+N23*L31+N25*L51 
R22=N22*L23+N23*L33+N25*L53 
R23=N22*L24+N23*L34+N25*L54 
R24=N22*L22+N23*L32+N25*L52 
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c 

WRITE (13,*) R11,R13,R22,R24 
K=(3*R11+R13+R22+3*R24) /8 
WRITE (11,2001) WN,K,CCRIT 
1 CONTINUE 
STOP 

1004 FORMAT (' ENTER NUMBER OF DATA') 

1006 FORMAT (' ENTER DELTAS AT. ' ) 

1007 FORMAT ('ENTER A3') 

2001 FORMAT (3E15.5) 

END 

C================================================= 

SUBROUTINE DSOMEG ( UK, WR, WI , OMEGA, CHECK) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION WR ( 6 ) , WI ( 4 ) 

CHECK=-1.0D+25 
DO 1 1=1,6 

IF (WR(I) .LT. CHECK) GO TO 1 
CHECK=WR ( I ) 

IJ=I 

1 CONTINUE 
OMEGA=DABS (WI ( IJ) ) 

IF (WI (IJ) .GT.O.DO) IJK=IJ 
IF (WI (IJ) .LT.O.DO) IJK=IJ-1 
RETURN 
END 

C================================================= 

SUBROUTINE NORMAL (T) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION T ( 6 , 6 ) , TNOR (6,6) 

CFAC=T(1, 1) **2+T(l,2 ) **2 
IF (DABS(CFAC) .LE. (l.D-10) ) STOP 4001 
TNOR (1, 1) =1 .DO 

TNOR (2, 1) =(T(1,1)*T(2,1)+T(2,2)*T(1,2)) /CFAC 
TNOR (3,1)=(T(1,1)*T(3,1)+T(3,2)*T(1,2)) /CFAC 
TNOR (4,1)=(T(1,1)*T(4,1)+T(4,2)*T(1,2)) /CFAC 
TNOR (5,1)=(T(1,1)*T(5,1)+T(5,2)*T(1,2)) /CFAC 
TNOR (6,1)=(T(1,1)*T(6,1)+T(6,2)*T(1,2)) /CFAC 
TNOR (1,2) =0.D0 

TNOR(2,2) =(T(2,2) *T(1,1) -T ( 2 , 1 ) *T ( 1 , 2 ) ) /CFAC 
TNOR(3,2) =(T(3,2) *T(1,1) -T (3 , 1 ) *T ( 1 , 2 ) ) /CFAC 
TNOR (4,2)=(T(4,2)*T(1,1)-T(4,1)*T(1,2)) /CFAC 
TNOR(5,2)=(T(5,2) *T(1,1) -T ( 5 , 1 ) *T ( 1 , 2 ) ) /CFAC 
TNOR (6,2)=(T(6,2)*T(1,1)-T(6,1)*T(1,2)) /CFAC 
DO 1 1=1,6 
DO 2 J=l,2 

T(I, J) =TNOR(I, J) 

2 CONTINUE 

1 CONTINUE 

RETURN 
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END 

C====================================================== 

SUBROUTINE MULT (TINV, A, T, A2 ) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION TINV (6,6),A(6,6),T(6,6),A1(6,6),A2(6,6) 
DO 1 1=1,6 
DO 2 J=l,6 
A1 (I, J) =0 .DO 
A2 (I, J) =0.D0 
2 CONTINUE 
1 CONTINUE 

DO 3 1=1,6 
DO 4 J=l, 6 
DO 5 K=l,6 

A1 (I, J) =A(I,K) *T(K, J) +A1 (I, J) 

CONTINUE 
CONTINUE 
CONTINUE 
DO 6 1=1,6 
DO 7 J=1 , 6 
DO 8 K=l,6 

A2 (I, J) =TINV(I, K) *A1 (K, J) +A2 (I, J) 

CONTINUE 
CONTINUE 
CONTINUE 
DO 11 1=1,6 
C WRITE (*,101) (A(I, J) , J=l, 6) 



11 


CONTINUE 
DO 12 1=1,6 




c 

12 


WRITE (*,101) 
CONTINUE 
DO 10 1=1,6 


(T(I, J) , J=l,6) 


C 

10 


WRITE (*,101) 
CONTINUE 


(A2 (I, J) , J=l, 6) 


C 


WRITE (*,101) 
RETURN 


A2 ( 1 , 1 ) 


101 


FORMAT (4E15 
END 


.5) 
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